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ABSTRACT 


A  computationally  efficient  scheme  for  multi  pie- source  location 
estimation,  based  on  the  Estimate-Maximize  (EM)  algorithm  [1],  is  presented. 
The  proposeo  scheme  is  optimal  in  the  sense  that  it  converges  iteratively  to 
tiie  exact  ^ldxi niurn  Livelihood  estimate  of  all  source  location  parameters 
simultaneously.  Versions  of  the  algorithm  that  incorporate  the  estimation  of 
tne  uiixnown  amplitude  attenuations  and  the  estimation  of  the  unknown  signal 
waveforms  are  also  presented. 


I.  INTRODUCTION 


Tiie  locdtion  of  a  radiating  point- source  can  be  determined  by 
Observation  of  its  signal  at  an  array  of  spatially  distributed  sensors.  The 
optimal ,  haxiniuni  Likelihood  (hL),  estimate  of  the  source  location  parameters 
(i.e.  bearing  ana  range)  is  achieved  by  maximizing  the  array  beamformer  output 
lZj  or,  alternatively,  by  cross-correlating  the  various  sensor  outputs  and 
obtain  the  ML  estimates  by  maximizing  a  weighted  sum  of  the  cross-correlation 
responses  witii  respect  to  a  pair  of  bearing  and  range  parameters  [3]. 

Tne  presence  of  several  signal  sources  drastically  complicates  the 
estimation  process.  To  obtain  the  ML  estimate  of  all  source  location 
parameters  jointly,  we  have  to  maximize  a  significantly  more  complicated 
fiiglTly  rion- 1  itiear  function  with  respect  to  K  pairs  of  unknown  location 
parameters,  where  K  is  the  assumed  number  of  signal  sources.  Of  course,  brute 
force  can  always  be  used  to  solve  the  problem,  evaluating  the  objective 
function  on  coarse  grid  to  roughly  locate  the  global  maximum,  and  then 
applying  the  Gauss  method  or  Newton- Raphson  or  some  other  iterative 
gradient- search  algorithm.  However,  when  applied  to  the  problem  in  hand, 
these  methods  tend  to  be  very  complex  and  computationally  time  consuming. 

Consequently,  approximations  to  the  ML  solution  and  various  ad  hoc 
sciiemes  for  multiple  source  localization  have  been  proposed  in  recent 
literature  (e.g.,  L4J  -  L20j);  however,  most  of  the  proposed  suboptimal 
localization  schemes  simplify  processor  structure  and  computations  at  the 
sacrifice  of  system  resolution  and  accuracy. 

Ill  this  report,  we  aevelop  an  iterative  scheme  for  multiple  source 
location  estimation  based  on  the  Estimate-Maximize  (EM)  algorithm.  However, 
uiiliXe  the  orute  force  gradient- search  iterative  algorithms,  the  proposed 
scheme  makes  an  essential  use  of  the  stochastic  system  under  consideration. 


The  heuristic  idea  is  to  decompose  the  sum  of  vector  signals  observed  at  the 
sensor  outputs  into  its  components,  and  then  apply  a  conventional  beamformer 
instrumentation  to  each  signal  component  to  obtain  the  bearing  and  range 
estimate  of  the  corresponding  source.  The  algorithm  iterates,  using  the 
resulting  parameter  estimates  to  better  decompose  the  observed  data  on  the 
next  iteration  cycle  and  thus  to  improve  the  next  parameter  estimates.  This 
computationally  attractive  algorithm  is  shown  to  be  optimal  in  the  sense  that 
it  converges  to  the  exact  KL  estimate  of  all  source  location  parameters 
siiiiul  taneously . 

Tne  organization  of  the  report  is  as  follows:  In  Section  II  we 
re-derive  the  EM  algorithm  following  the  considerations  in  [1],  and  then  we 
specialize  to  the  Linear-Gaussian  case.  In  Section  III  we  apply  the  algorithm 
to  the  multiple  source  location  estimation. 


11.  MAXIMUM  LIKELIHOOD  ESTIMATION  AND  THE  ESTIMATE -MAX I MIZE  ALGORITHM 
Let  T  be  a  vector  random  variables  possessing  the  probability 
density  (zf  ,  where  Q  is  an  open  subset  of  the  K-dimensional 

Euclidean  space.  The  ML  estimate  of  9  ,  given  an  observed  W,  is 
outained  by  maximizing  the  log-likelihood  function 


9  6 


(1) 


If  the  vector  9  contains  several  unknowns,  such  as  in  the  multiple 
source  location  problem,  and  since  f  Y'  ^ generally  a  highly 

nun- linear  function  of  9  ,  the  maximization  required  in  (1)  tends  to  be  very 
cumplex. 


•  3  • 


Suppose  the  data  vector  Y  can  be  viewed  as  being  "Incomplete*,  and  we 
can  specify  some  "complete*  data  X  that  Is  related  to  Y  by 


=  Y 


where  j-l/'-Jis  a  non-invertable  (many-to-one)  transformation.  In  the  multiple 
source  location  problem,  the  "complete"  data  Xcould  be  the  observation  of  the 
various  source  signals  separately,  where  the  "incomplete"  (observed)  data  Y  is 
ti»e  s^i  of  the  signal  contributions  from  the  various  sources.  As  pointed  out 
before,  the  Y  model  may  be  complicated  to  work  on  directly,  in  which  case 
reference  to  theX^odel  might  be  very  useful. 

For  all  "complete"  data  realizations  X  such  that 


(3) 


wnere  /^x^f^is  the  probability  density  of  X  ,  and  fx/Y:  ^ 
conditional  probability  density  of  X  given  that  Y?  ^  .  Taking  the  logarithm 
on  both  sides  of  (3),  we  obtain 


TdKing  tne  conditional  expectation  given  Y:  ^  for  a  parameter  value 
(i.e.,  multiplying  both  sides  of  (4)  integrating  with 
respect  to  X  over  {  obtain 


^^3  fy  =  f  /”  h  ix  -'A-‘  ^  j 
-flh  /x/J'  9 
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wliere  we  note  that  for  a  fixed  ^  ,  the  left  hand  side  of  (5)  is  a  constant  and 
hence  it  is  unaffected  by  the  conditional  expectation  operation.  Define  for 
convenience 

E  (  /yfij  ^  e‘]  (6) 

P  ^  ff  ?']  (7) 

With  tnese  definitions,  (5)  reads 

^^5  £?0  (8) 

Using  the  Jensen's  inequality,  it  can  easily  be  verified  that 
P f)  ^  '  Hence 

(pPpP'^>  ==^  <9) 

Tne  relation  in  (y)  fonn  the  basis  to  the  EM  algoritin.  The  algon■t^^l 
starts  with  an  initial  guess  and  let  be  defined  inductively  by 


b 


9 


(n^i) 


77^0,1,2,... 


(10) 


Tne  Jensen's  inequality  asserts  that  for  any  two  probability 
densities and ^Oi )defined  over  fX 


^{(i)  fcj 

II 


-fl 
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Since  the  value  of  B  that  maximizes  ,  then 

accorainy  to  (9)  each  iteration  of  the  algoritJm  increases  the  likelihood. 
Hence,  under  the  usual  regularity  conditions,  the  algorithm  converges  to  the 
maxiiiiun.  likelihood  estimate,  i.e.  .  The  rate  of  convergence  of 

the  algorithm  is  exponential  (see  [1]),  depending  on  the  fraction  of  the 
covariance  of  the  “complete"  data  that  can  be  predicted  using  the  observed 
("incomplete")  data.  If  that  fraction  is  small,  the  rate  of  convergence  tends 
to  be  slow,  in  which  case  one  could  use  standard  numerical  methods  to 
accelerate  the  algorithm. 

It  is  important  to  observe  that  the  EM  algoritln  is  not  uniquely 
defined.  The  transformation  relating  X  to  V  can  be  any 

noi.-invertable  transformation.  Obviously,  there  are  many  possible  "complete" 
data  specifications  that  will  generate  the  observed  data.  Thus,  the  EK 
algorithm  can  be  implemented  in  many  possible  ways.  The  final  outcome,  which 
is  the  ML  estimate,  is  completely  unaffected  by  the  way  in  which  H  is 
specified  (i.e.  the  choice  of  "complete"  data);  however,  the  choice  of  Hmay 
critically  affect  the  complexity  and  rate  of  convergence  of  the  algorithm,  and 
unfortunate  choice  of  H  may  yield  a  completely  useless  algoritim.  Hence, 
given  a  parameter  estimation  model,  the  practically  important  question  that  is 
left  open  is  how  to  find  the  computationally  most  efficient  implementation  of 
the  algorithm. 


The  Linear-Gaussian  Case 

Suppose  that  HX’Y »  H  matrix  and  ^possesses 

the  following  multivariate  Gaussian  probability  density 


1 


(11) 


•  6  ^ 

wliere  if  V  is  real-valued,  2  If  V  is  complex- valued,  and  the 
superscript  P-  denotes  the  conjugate-transpose  operator.  We  shall  refer  to 
this  case  as  the  Linear-Gaussian  case.  Taking  the  logarithm  of  (11),  we 
obtain 

'  “  ^  ^  >3^ /"f)  (12) 

wiiere  trL  J  stands  for  the  trace  of  the  bracketed  matrix.  Substituting  (12) 
into  (6) ,  we  obtain 

wr«r.  /  =  f  i  '  }  .and  /x' rf ^  j  f '  }  . 

Using  well-known  results  from  linear  estimation  ([23],  Chapter  4),  we  obtain 

Vs  02^/f') H->^  (14) 

fi*'-  [I- rc^n-HlR/i')'^ 

wiiereP/flis  the  "Kalman  Gain*  defined  by 

r/'e)z 


(16) 
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Substituting  (14)  and  (15)  into  (13),  the  function  ^<V^/Orequired  for 
the  EK  algorithm  is  given  in  a  closed  form.  We  note  that 
and  Pc’j the  sau)e  dependence  on  f  .  Maximiring  t')  with 
respect  to  t  is  the  same  as  maximizing  i)  with  respect  to  e  . 

Hence,  tiie  EK  algoritlim  requires  the  ML  solution  in  the  X model  which  might  be 
sigtiificantly  easier  to  obtain  than  the  direct  hL  solution  in  the  V  model . 

If  ^  constant  matrix,  f)  assumes  the  simplified  form 

^  z  ^  y  i  [/*  ) /e  W/pJ  (17) 

wnere  accounts  for  all  terms  that  are  independent  of  b.  Substituting 
(17)  ifito  (10),  the  Estimate  (conditional  expection)  step  and  Maximize  step  of 
the  algorixhn  are  given  by: 

E-step:  Compute 


(16) 


M- step: 

f 

III.  ARRAY  PROCESSING  VIA  THE  EK  ALGORITHM 

The  basic  system  of  interest  consists  of  K  spatially  distributed  signal 
sources  and  an  array  of  M  spatially  distributed  sensors  as  illustrated  in 
Figure  1.  Assuming  perfect  propogation  conditions  in  the  medium  between 
sources  and  receivers,  the  actual  waveform  observed  at  the  m^^  sensor  output 
is  given  by 

1C 

»*/ 


(19) 


(20) 


where  the  source  signal,  'y7^li)\i  the  additive  noise,  T/^yy,  is 

the  travel  time  of  the  signal  wavefront  from  the  source  to  the 
sensor,  and  o(.  are  the  amplitude  attenuations. 

Infonaation  concerning  the  various  source  location  parameters  can  be 
extracteo  b>  measuring  the  various  .  In  the  passive  case,  one  can  only 
measure  the  travel  time  differences,  obtainable  by  selecting  one  sensor  as  a 
reference  ano  comparing  its  output  with  that  of  every  other  sensor.  If  we  let 
sense*  Af  to  be  the  reference  and  set  measures  the  travel 

time  difference  of  the  signal  wavefront  to  the  [yy^,M  )  sensor  pair. 

To  simplify  the  exposition,  suppose  that  the  various  signal  sources  are 
relatively  far-field  so  that  the  observed  signal  wavefronts  are  essentially 
planar  across  the  array.  If  we  further  suppose  that  the  array  sensors  are 
co-linear,  then 

^  'fi  (21) 

where  is  the  spacing  between  sensors and  <  is  the  velocity  of 
propogation  in  the  medium,  and  is  the  angle  of  incident  of  the  k^^  signal 
wavefront  with  respect  to  the  array  baseline. 

In  this  setting,  the  estimation  problem  can  be  stated  as  follows:  Given 
the  ooserveo  data  .  ^ihd  the  H  estimate  of  1^^  ...  . 

We  shall  find  it  convenient  to  work  with  the  parameters  9/^  •  Ctfi 
Since  M.  estimation  commutes  over  non-linear  transformation,  we  can  first 
estimate  the  9/^  and  then  translate  to  the  . 

Assuming  that  the  5  are  perfectly  known  to  the  observer,  and 

that  the  are  realizations  from  uncorrelated  zero-mean  spectrally 

white  Gaussian  processes,  the  log- likelihood  function  is  given  by 

/W  y  ^ 

^03  jr  'I.  J’’ 

T 

*  4 


where  =  dm /<  ,  Is  the  spectral  level  of  ,  [T*^T>  3  Is  the 
observation  interval,  and  ^  is  a  normalizing  constant.  The  result  in  (22)  is 
a  straightforward  multi-channel  extension  of  the  known  (deterministic)  signal 
in  white  Gaussian  noise  problem  ([24],  Chapter  4).  Thus,  the  M.  estimate  is 
the  solution  to  the  following  problem; 


M 


Tc 


A 


(23) 


Ignoring  terms  that  are  independent  of  f  ,  (23)  reduces  to 

T/ 


J 

(  ry:  1 

/ 

n 

-  HI  pi  =^i 

i-l  Ol 


7". 


In  many  situations  of  practical  interest,  the  number  of  signal  sources 
is  also  an  unknown  parameter  to  be  estimated,  in  which  case  several  criteria 
to  determine  IC',  based  on  the  ML  estimate  of  the  '6  ,  are  presented  in 
(L20J  -  LZcj).  Thus,  to  obtain  the  Ml  estimate  of  the  jointly  with  the 

estimate  of  the  number  of  signal  sources,  we  have  to  solve  a  non-linear 
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optimization  problem  in  K  unknowns,  for  which  a  closed  form  analytical 
solution  cannot  be  found. 

Wb  further  note  that  the  optimization  problem  stated  above  assumes  prior 
knowledge  of  the  and  practice,  this  Is  apt  to  be 

unrealistic.  One  is  unliiiely  to  have  exact  prior  knowledge  of  the  amplitude 
attenuations  ana  detailed  description  of  the  signal  waveshapes  may  be 
similarly  incomplete,  in  which  case  to  obtain  the  H  estimate  of  the 
we  must  maximize  the  expression  in  (24)  with  respect  to  all  the  unknowns  in 
the  problem.  The  effect  of  unknown  attenuation  factors  can  be  eliminated  from 
the  likelihood  equation  by  observing  that  for  pre-specified  ^ ,  (23) 
becomes  a  weighted  linear  least  squares  problem  In  the  ,  for  which  there 

is  a  closed  foroi  solution.  Thus,  we  can  substitute  the  ^  by  their 

weigiited  least  squares  estimates  and  obtain  a  somewhat  more  complicated 
functional  to  be  maximized  with  respect  to  the  remaining  unknowns.  However, 
the  effect  of  unknown  signal  waveshapes  cannot  be  eliminated  that  easily  and 
hence,  the  required  maximation  becomes  exceedingly  more  complicated. 

Having  the  EH  method  in  mind,  we  would  like  to  simplify  the  maximization 
associated  with  the  direct  ML  approach.  To  apply  the  algorithm  to  the  problem 
in  hand,  the  first  step  is  to  specify  the  "complete"  data.  A  natural  choice 
of  the  "complete"  data  is  given  by  decomposing  the  observed  signals 

'Xu  ‘  t  H  )  (25) 

Where  the  chosen  to  be  realizations  from  uncorrelated  zero-mean 

white  Gaussian  processes  with  spectral  levels  of  require 

that 

ic 

7  ^4^  /f  ^ 


(26) 
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then 


k: 


I  yt 


*r) 


U) 


(t> 


(27) 


Concatenating  the  M  equations  In  a  vector  form,  we  obtain 

kz!  " 


where 


(29) 

and 

^j(i)  ^  \(i)  (30) 


Equation  (28)  can  be  rewritten  in  the  form 


H  yfi)-  jH) 


(31) 
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nthe  rfe 


y(i)  * 


f  ) 


(32) 


onu 


H= 

■ - > - ' 


(33) 


where  I  is  the  MXM  identity  matrix.  We  denote  //^)  the  “complete"  data, 
and  '^(()the  “incomplete"  (observed)  data.  Note  that  /(Ois  composed  of  the  K 
components  y/; /i)  ^  ,  where  V^/t)  is  the  vector  signals  received 

from  tie  source  alone. 

tonsioer,  for  tlie  moment,  M.  estimation  using  the  "complete"  data. 

Since  the  components  of  //V)are  uncorrelated  Gaussian  random  processes,  the 
loy-lixelihooQ  of is  the  sum  over  i  and  w  of  the  log- likelihood  of 
the  ,/<)  which,  in  turn,  is  given  by 


a  - 


Suppose  that  the  are  known  a-priori  so  that  we  "only"  have  to 

r 

estimate  the  K  sets  ^  _ maximization  of  (24) 

with  respect  to  all  unknowns  is  a  rather  complex  maximization  problem.  The 
maximization  of  (34),  however,  is  a  much  simpler  problem.  This  is  since  the 
expression  in  (34)  can  be  decomposed  into 


(3A) 


(35) 
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where 


^j(  2j 

^  Iv  i  ^  77  ^  fi'  fr^  )'\  ^ 

£>  Ia^  ^  ~  tm  ^ 

Ttie  expression  in  (36)  depends  only  on  ^  ^ii  >  ••-  *<1^  • 
Hence,  the  maximization  of  (34)  with  respect  to  the  K-set  unknowns  can  be 
decomposea  into  the  K  separate  maximizations 


(36) 


/  - 
May  L. 


A/l 


f  r  ( 

i  r. 


T. 


(37) 


The  maxtmization  in  (37)  is,  in  fact,  the  maximization  problem 
associoteo  with  the  M.  estimation  of  the  bearing  parameter  of  a  single  source 
in  the  presence  of  unknown  amplitude  attenuations. 

The  EK  algorithm  is  directed  at  finding  a  value  of  the  parameters  that 
maximize  (24),  however,  it  does  so  by  making  an  essential  use  of  the  solution 
to  (37). 

Since  the  “complete"  data  is  Gaussian,  and  the  transformation 
relating  Y  /i')  to  y  /t  1  is  linear,  we  can  use  the  version  of  the 
algorithm  developed  for  the  Linear-Gaussian  case.  Thus,  the  M-step  of  the 
algorithm  is  given  by  (34)  (since  fO  and  ^  ) 

assumes  the  same  form),  where  the  components  of  ///)are  substituted  by 

f  [ yf*'>/HyN):  f '"J 


(38) 
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Tht  conditional  expectation  required  in  (38)  can  easily  be  obtained 
usiitj^  (lb).  The  resulting  algorit^  is  given  by: 

E-step 


Coiiipute 


f-! 


(39) 


M- step 


For 


h!,  - 


T/ 


,  o(j.^  , Ti 


(40) 


Perhaps  the  must  stricking  feature  of  the  algoritim  has  already  been 
indicateo  before.  The  algorithm  decouples  the  complex  multiple-parameter 
maximization  associated  with  the  direct  FL  procedure  into  K  separate  M. 
maximizations  as  illustrated  in  Figure  2.  The  extension  to  bearing  and  range 
estimotion  is  straightforward.  The  basic  scheme  is  still  illustrated  by 
Figure  2,  where  now  each  ML  processor  requires  the  maximization  with  respect 
to  a  pair  of  bearing  and  range  parameters.  Thus,  the  complexity  of  the 
algorithm  is  completely  unaffected  by  the  number  of  signal  sources.  As  the 
number  of  sources  increases,  we  have  to  increase  the  number  of  FL  processors 
in  parallel;  however,  each  processor  is  maximized  separately. 
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Since  the  algorithm  is  based  on  the  EM  method,  it  must  converge  to  the 
exact  ML  estibtate  of  all  source  location  parameters  simultaneously. 

We  note,  in  passing,  that  according  to  (26),  the  must  satisfy  the 
constraint 

K 

I  -  /  (41) 

/--/ 


but  otherwise  they  are  free  variables.  The  choice  of  does  not  affect  the 
final  estimate  (at  the  point  of  convergence);  however,  it  can  be  used  to 
control  the  rate  of  convergence  of  the  algoritlm. 

We  now  consider  the  M-step  of  the  algorithm  in  more  details.  If  we  set 
the  Derivative  of  the  expression  in  (40)  equals  to  zero,  we  obtain 


(42) 


Since  the  second  derivatives  with  respect  to  the  ®  negative 

oefinite  matrix,  we  have  therefore  expressed  the  optimal  choice  of  the  as 
a  function  of  .  Substituting  (42)  into  (40)  and  following  straightforward 
algeora  manipulations,  the  M-step  of  the  algorithm  reduces  to 


uh- 


(A3) 


Wl  :/ 


T, 


1  « 

/ r/vif-  tz.  ) *f i 

T. 


(44) 
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r. 


Ill  most  situations  of  practical  interest,  we  can  assume  that 
J  is  independent  of  ,  in  which  case  the  above  pair 

of  equations  reduces  to 

^  j  ^  .  T> 


T I 

r 

r-y  i 


(45) 


0^ 


'.y*^ 


^  J  si' 

T  i 


(46) 


Tf 


The  tenn  f  ft)  can  be  generated  by  passing 

through  a  filter  matched  to  fi).  Thus,  is  obtained  by  maximizing 

a  squareo  ana  weighted  sum  of  a  bank  match  filter  outputs. 


Modified  EH  Algorittiii 

The  E>.  theory  allows  us  to  substitute  the  f^-step  of  the  algorithm 
(Equation  (40))  by  the  following  two-step  maximization: 

T. 


y 


fa  7  ' 


-  j  ==#> 


(47) 


r, 


T. 


5 


T; 


(48) 
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Assuming,  as  before,  that 


and  (4b)  reouce  to 


T. 


Tf 

Is  independent  of  ,  (47) 


z  K'-Vo^  fi-r^  eiUt  — ^ 

m-/  /A  )  ^  k 


^  7. 


y 


(y"’  l) 


7/ 


(49) 


(50) 


We  note  that  the  solution  to  (49)  is  obtained  by  maximizing  a  weighted 
linear  sum  of  a  bank  match  filter  outputs. 

Tne  bioximization  of  (47)  followed  by  the  maximization  of  (48)  is 
generally  not  equivalent  to  the  maximization  required  in  (40).  However,  since 
the  condition  in  (9)  is  still  satisfied,  each  iteration  increases  tite 
likelihooa,  and  the  modified  algorithm  converges  to  the  desired  results.  By 
replacing  (4b)  by  (49)  ,  we  therefore  simplify  processor 

structure  ano  computations  in  the  expense  of  possibly  very  moderate  decrease 
in  the  rate  of  convergence  of  the  algoritlm. 

Another  alternative  is  to  maximize  the  expression  in  (40)  first  with 
respect  to  the  (using  fcj'"’)  and  then  with  respect  to  9^^  .  The  H-step  of 
the  algorithtTi,  in  that  case,  is  given  by 

-7Tf  ; 

T/ 


(51) 
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f  f- 

I  -77 —  y/^Vi)  J/ /V-C.  an 

r.  »•■'  J  t-  ■’^  t 

‘'i  T. 

We  note  that  (49)  followed  by  (50)  is  not  the  same  as  (51)  followed  by 
(bii).  However,  for  the  same  reason,  both  algorithms  converge  to  the  desired 
ML  solution. 

Unknown  Signal  haveforms 

The  algorithm  and  its  modified  version  developed  in  the  previous 
sections  critically  depends  on  exact  prior  knowledge  of  the  5^ (t).  In 
practice,  one  is  unlikely  to  have  detailed  prior  knowledge  of  the  signal 
waveforms  ,  in  which  case  they  must  be  estimated  jointly  with  the  and 

tny  •  Following  the  same  considerations  as  in  the  known  signal  case, 

the  resulting  algorithm  is  given  by 
E-step 

Compute 


►t-  step 


For  /  =  6 


T, 


9k ^  U(i)^  o/k,j...  9^^^ 

-  p/vy-  ft  )/il  — spi)^  v/;"’ 

T. 


(53) 


(54) 
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Again,  the  EM  algorithm  converts  the  extremely  complicated  maximization 
problem  associated  with  the  direct  M.  procedure  Into  K  separate  maximizations 
as  suggested  by  Figure  2.  Each  maximization  Is,  in  fact,  the  maximization 
problem  associated  with  the  M.  estimation  of  the  bearing  parameter  of  a  single 
source  of  unknown  waveshape  in  the  presence  of  unknown  amplitude  attenuations. 
Ignoring  end>effects,  the  expression  in  (54)  to  be  maximized  can  be 

written  in  the  form 

T( 

if' I  H  (ss) 

/I  w*/  ‘‘ 

T, 


To  find  a  function  Si(t)  that  maximizes  the  integral  given  in  (55),  it 
is  sufficient  to  find  .^^^ithat  maximizes  the  integrand.  If  we  set  the 
derivative  of  the  integrand  with  respect  to  SkU)  equals  to  zero,  we  obtain 
after  some  obvious  manipulations 

A) 

_  (56) 


U) 


Since  the  second  derivative  of  the  integrand  with  respect  to  is 
negative,  we  have  found  the  optimal  choice  of  expressed  in  terms  of  the 

remaining  unknowns  in  the  problem.  Substituting  (56)  into  (54),  the  required 
maximization  can  be  carried  out  as  following: 
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tAc^y  )  ^  7 


n<-.  /  (i!  ^/i£ 


Ailj  ... 


'f 

T.  / 


■»  S’/  ,  •^k,  ,  ...  Uf. 


y> 


r 


(57) 


(58) 


yyi-/ 


We  note  that  at  convergence,  equation  (56)  yields  the  ML  estimate  of 
tne  as  well.  This  information  is  of  considerable  interest  in  many 

applications. 

The  maximization  required  in  (57)  is  still  rather  complicated.  However, 
the  modified  EM  approach  can  be  used  to  further  simplify  the  computations  by 
maximizing  (&4)  first  with  respect  to  (P^  ,  then  with  respect  to  and 

finally  with  respect  to  the  .  The  M-step  of  the  algorithm,  in  that  case, 
taxes  tne  form: 


lA^y 

&k 


Al 

7  J 


*».*/  (»! 


^ke 


yk^'’n^y^ei)x,';  Vi’ ft  fiVi 


T, 


fyii') 

h 


2 

z  {42Y/fY^^ 

7»^s/ 


(60) 
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(59) 


(61) 
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In  this  setting  we  only  have  to  perform  a  one-parameter  maximization  at 
each  iteration  cycle.  The  expression  in  (59),  to  be  maximized,  can  be 
rewritten  in  the  form: 


(62) 


wnere  Ju;')  is  the  Fourier  transform  of  and  VI' is  the  signal 

frequency  band.  The  expression  in  (62)  is  the  array  beamformer,  implemented 
in  eititer  tne  time  or  the  frequency  ooroain.  Thus,  the  M-step  of  the  algorit)ri 
essentially  consists  of  maximizing  K  beamformers  in  parallel  as  illustrated  in 
Fi gure  o . 


Alternatives  to  Ute  K-step  can  be  obtained  by  changing  the  order  of 
maximizations.  For  example,  we  can  apply  the  various  algorithms  derived  for 
tlie  xnuwn  signal  case,  where  1/ /t)  is  substituted  by  its  current 
estimate  and  then  use  (60)  to  update  the  estimate  of  -f/t  )  .  The 

various  olgoritf«is  may  have  slightly  different  convergence  properties; 
however,  all  of  which  converge  to  the  exact  ML  estimate  of  all  source  location 
parameters  simultaneously. 
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COmLHlS 

1)  Tfie  extension  of  the  algorithms  to  bearing  and  range  estimation  and 
arbitrary  (out  xnown)  array  configurations  is  straightforward.  2)  The 
derivation  of  tne  algoritlims  for  the  case  in  which  the  M  5  are  modelled 
as  sample  functions  from  uncorrelated  zero-mean  Gaussian  random  processes  is 
presented  in  [25],  d)  In  the  assumed  model,  we  have  ignored  the  phase  shifts 
caused  by  scattering  and  reflection  phenomena.  These  effects  can  be  taken 
into  account  by  considering  the  in  (20)  to  be  the  complex  envelope  of 

the  receiveo  signals,  in  which  case  the  are  complex- valued  (magnitude 

and  phase)  amplitude  attenuations.  The  extension  of  the  algorithms  to  the 
complex  case  is  straightforward  (note  that  the  Linear-Gaussian  case  has  been 
developed  for  complex  processes  as  well).  4)  The  proposed  scheme  can  be 
extended  to  time  delay  and  location  estimation  in  multipath  environment. 


23 


ACKNOWLEDGEMENTS 


This  stucly  has  been  supported  by  The  Naval  Air  Systems  Command  under 
contract  Nuu014-8b-K-0272.  The  authors  wish  to  thank  Ms.  Cindy  Leonard  for 
her  excellent  secretarial  assistance. 


-  24  - 


REFEKENCES 


Lij  A.F.  Deupster,  K.M.  Laird,  and  D.B.  Rubin  (1977):  "Maxiinuni  Likelihood 
from  Incomplete  Data  via  the  EM  Algorithm,"  Ann.  of  the  Royal  Stat.  Soc.. 
pp.  l-3b. 

IZj  W.J.  bangs  and  P.M.  Schultheiss  (1973):  "Space-Time  Processing  for 
Optimal  Parameter  Estimation,"  in  Signal  Processing,  J.W.R.  Griffiths, 
P.L.  Stocklin,  and  C.  Van  Schooneveld,  Eds.  New  York:  Academic,  pp. 
577-590. 

L3j  K.k.  Hahn  (1975):  "Optimum  Signal  Processing  for  Passive  Sonar  Range  and 
Bearing  Estimation,"  J.  Acoust.  Soc.  Amer.,  Vol.  58,  pp.  201-207. 

L4j  Borgiotti,  G.V.  and  Kaplan,  L.J.  (1979):  "Superresolution  of 

Uncorrelated  Interference  Sources  by  Using  Adaptive  Array  Techniques," 
IEEE  Trans,  on  Antennas  and  Propagation,  Vol.  27,  pp.  842-845. 

Lbj  Cantoni,  A.  and  Godara,  L.C.  (1980):  "Resolving  the  Directions  of  the 
Sources  in  a  Correlated  Field  Incidenting  on  an  Array,"  Journal  of 
Acoustical  Society  of  America,  Vol.  67,  pp.  1247-1255. 

Loj  oe  Figueiredo,  k.J.P,  and  Gerber,  A  (1983):  "Seperation  of  Superimposed 
Signals  Dy  Cross-Correlation  Method,"  IEEE  Trans,  on  Acoustics,  Speech 
and  Signal  Processing,  Vol.  31,  pp.  1084-1089. 

l7j  Nehorai,  A.,  Su,  G.  and  Morf,  M.  (1983):  "Estimation  of  Time  Difference 
of  Arrival  for  Multiple  ARHA  Sources  by  Pole  Decomposition,"  IEEE  Trans, 
on  Acoustics,  Speech  and  Signal  Processing,  Vol.  31,  pp.  1478-1491. 

Loj  Paulraj,  A.  and  Kailath,  T.  (1985):  "Direction  of  Arrival  Estimation  by 
Eigen  structure  With  Unknown  Sensor  gain  and  Phase,"  Proc.  IEEE  Int.  Conf. 
on  Acoustics,  Speech  and  Signal  Processing,  (Tampa,  FL). 

L9j  Porat,  B.  and  Friedlander,  B.  (1983):  "Estimation  of  Spatial  and 
Spectral  Parameters  of  Multiple  Sources,"  IEEE  Trans,  on  Information 
Theory,  Vol.  29,  pp.  412-425. 

LlOj  Reodi,  S.S.  (1979):  "Multiple  Source  Location  -  A  digital  Approach," 

IEEE  Trans,  on  Arespace  and  Electronic  Systems,  Vol,  15,  pp.  95-105. 

Lilj  Schmidt,  R.O.  (1981):  "A  Signal  Subspace  Approach  to  Multiple  Emitter 
Location  and  Spectral  Estimation,"  Ph.D.  dissertation,  Stanford 
University,  CA. 

L12J  Sclimidt,  R.O.  and  Frank,  R.E.  (1983):  "Multiple  Source  DF  Signal 
Processing:  An  Experimental  System,"  Proc.  Symposium  on  Antennas 
Applications,  (Montecello,  IL). 

Ll3j  Scnweppe,  F.C.  (196B):  "Sensor  Array  Data  Processing  for  Multiple  Signal 
Sources,"  IEEE  Trans,  on  Information  Theory,  Vol.  14,  pp.  294-305. 


-  25  - 


[14j  Su,  G.  and  Morf,  M.  (1962):  "The  Signal  Subspace  Approach  for  Multiple 
Emitter  Location,"  Proc.  16th  Asilomar  Conf.  on  Circ.  Syst.  and 
Conputers,  (Pacific  Grove,  CA),  pp.  336-340. 

Llbj  Su,  G.  and  Morf,  M.  (1983):  "The  Signal  Subspace  Approach  for  Multiple 
Wideband  Emitter  Location,"  IEEE  Trans,  on  Acoustics,  Speech  and  Signal 
Processing,  Vol.  31,  pp.  1503-1522. 

Lloj  Wang,  H.  and  Kaveh,  M.  (1984):  "Estimation  of  Angles- of-Arrival  for 

Wideband  Sources,"  Proc.  IEEE  Int.  Conf.  on  Acoustics,  Speech  and  Signal 
Processing  (San  Uiego,  CA),  pp.  7.5 .1-7 .5.3. 

Ll7j  Wax,  h.  and  Kailath,  T.  (1983a):  "Optimum  Localization  of  Multiple 

Sources  by  Passive  Arrays,"  IEEE  Trans,  on  Acoustics,  Speech  and  Signal 
Processing,  Vol.  31,  pp.  1210-1218. 

LlBJ  Wax,  M.  (1985):  "Detection  and  Estimation  of  Superimposed  Signals," 
Ph.D.  dissertation,  Stanford  University,  CA. 

Ll9j  Wax,  M.  and  Kailath,  T.  (1984b):  "Decentralized  Processing  in  Passive 
Arrays,"  Proc.  IEEE  Int.  Conf.  on  Acoustics,  Speech  and  Signal 
Processing,  (San  Diego,  CA),  pp.  40.7.1-40.7.4. 

l20j  Wax,  K. ,  Shan,  T-J.  and  Kailath,  T.  (1982):  "Location  and  Spectral 
Density  Estimation  of  Multiple  Sources,"  Proc.  16th  Asilomar  Conf.  on 
Circ.  Syst,  and  Computers,  (Pacific  Grove,  Ca),  pp.  322-326. 

L21j  Wax,  M.  and  Kailath,  T.  (1983b):  "Estimating  the  Number  of  Signals  by 
Information  Theoretic  Criteria,"  ASSP  Spectral  Estimation  Workshop  II. 
(rampa,  FL),  pp.  192-196. 

L22j  Wax,  M,  and  Kailath,  T.  (19b4a):  "Determining  the  Number  of  Signals  by 
Information  Theoretic  Criteria,"  Proc.  IEEE  Int.  Conf.  on  Acoustics, 
Speecn  and  Signal  Processing,  (San  Diego,  CA),  pp.  6. 3. 1-6. 3. 4. 

t23j  Gelo,  A.  (1974):  Applied  Optimal  Estimation,  Cambridge:  M.I.T.  Press. 

L24j  H.L.  Van  Trees  (1968):  Detection,  Estimation  and  Modulation  Theory, 
Part  I,  New  York:  Wiley. 

L2t)j  Feder,  H.  and  Weinstein,  E.  (1985):  "Optimal  Multiple  Source  Location 
Estimation  Via  the  EM  Algorithm,"  Proc.  IEEE  Int.  Conf.  on  Acoustics, 
Speech  and  Signal  Processing,  (Tampa,  FL)  pp.  1762-1765. 


-  26  - 


Figure  1 
Figure  Z 
Figure  3 


FIGURE  CAPTIONS 
Array- Source  Geometry 

Multiple  Source  Localization  via  the  EM  Algoritlm 

Localization  of  Multiple  Sources  Having  Unknown 
Waveforms 
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